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Abstract  — A  computer  model  was  designed  based  on  the 
concept  of  common  drive  which  suggests  that  motor  units 
(group  of  muscle  fibers  and  the  single  alpha-motoneuron  that 
innervates  them)  in  a  muscle  are  controlled  by  a  common 
input  to  the  entire  motoneuron  pool.  Where  possible,  the 
model  utilized  experimentally  determined  data  and 
supplemented  these  with  findings  reported  in  the  literature.  It 
was  validated  by  matching  the  simulated  mean  firing  rates, 
power  spectra,  and  compound  muscle  force  outputs  to  that 
produced  by  data  from  the  Tibialis  Anterior  muscle.  The 
model  was  implemented  using  Matlab’s  ®  SIMULINK  ® 
tool.  In  this  form,  the  model  allows  easy  modification  of 
parameters  to  allow  for  virtual  experimentation  that  would 
otherwise  be  impossible  with  human  or  animal  models.  The 
developed  model  was  used  to  evaluate  a  commonly  used 
technique,  spike-triggered  averaging  (STA),  to  estimate  the 
twitch  force  of  an  individual  motor  unit.  It  was  concluded  that 
STA  has  the  potential  to  produce  valid  estimates  only  at  firing 
rates  below  3  pulses  per  second  which  are  physiologically 
unfeasible.  Simulations  suggest  that  the  effects  of  common 
drive  on  reliable  MU  twitch  estimation  may  not  be  as 
extensive  as  initially  expected.  Additionally,  hypotheses 
regarding  the  effect  of  various  mechanical  characteristics 
under  certain  physiological  paradigms  such  as  hand 
dominance  or  fatigue  on  the  electrical  properties  can  be 
investigated  using  the  model. 
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I.  INTRODUCTION 

Computer  models  provide  a  means  to  investigate 
phyisological  systems  which  are  otherwise  forbidding  to 
experiment  with  either  due  to  the  difficulty  in  accessing  or 
modifying  defining  parameters  of  the  system.  The  present  model 
was  based  on  the  model  of  Fuglevand  et  al.  [7]  which  simulated 
isometric  force  from  a  model  that  predicted  recruitment  and  firing 
times  in  a  pool  of  motor  units.  It  differed  on  a  few  points:  The 
present  model  explicitly  assumed  a  design  that  reflected  the 
concept  of  common  drive  in  that  each  motor  unit  received  a 
common  input-  the  common  drive-  and  an  input  that  was  unique 
to  that  unit  and  not  shared  by  others-the  noise  signal  for  the 
purposes  of  this  model.  The  current  model  also  differed  from  the 
Fuglevand  model  in  various  firing  rate  characteristics  which  were 
designed  to  reflect  our  observations  in  various  but  most 
specifically  the  Tibialis  Anterior  muscle  [6]:  the  initial  firing  rate 
of  a  motor  unit  was  dependent  on  its  recruitment  threshold;  the 
excitation-firing  rate  profiles  were  not  composed  of  a  linear  region 
followed  by  a  saturation  region,  but  were  piecewise  linear 
throughout  the  excitation  range;  the  excitation-firing  rate  profile 
of  each  motor  unit  had  a  different  shape.  Furthermore,  the  pulse 
trains  were  estimated  from  the  firing  rate  estimates  by  Integral 
Pulse  Frequency  Modulation  rather  than  solving  for  interfiring 
intervals  and  hence  firing  times.  The  variability  in  interfiring 
intervals  was  achieved  by  the  inclusion  of  the  noise  signal  at  the 
input  rather  than  addition  of  noise  to  the  interfiring  intervals. 
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Figure  1  Model  of  force  production  in  the  muscle.  Based  on  the  concept  of  common  drive,  each  motor  unit 
receives  the  common  input  in  addition  to  a  noise  signal.  The  drive  is  translated  into  a  firing  rate  based  on  the 
characteristics  of  the  motor  unit,  which  are  tightly  associated  with  its  rank.  Integral  Pulse  Frequency  Modulation 
stage  generates  the  firing  instances.  The  impulse  response  to  the  firings  is  the  twitch  waveform  defining  the  force 
contribution  of  the  motor  unit.  Muscle  force  is  the  summation  of  contributions  of  all  motor  units. 
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II.  THE  MODEL 


C.  Firing  Train  Generation 


A.  Underlying  Theory 

The  observation  that  the  firing  rates  of  motor  units  fluctuate 
in  unison  with  essentially  no  time  delay  between  them  has  lead  to 
the  concept  of  common  drive  [4,  5],  This  finding  suggested  that 
the  CNS  has  evolved  a  relatively  simple  strategy  for  controlling 
motor  units.  Rather  than  controlling  the  activity  of  each  motor 
unit  separately,  the  CNS  appears  to  control  the  excitation  to  the 
motomeuron  pool.  The  common  drive  received  by  all  the  motor 
units  in  the  pool  gets  translated  into  individual  firing  patterns  for 
each  motomeuron  by  the  input/output  characteristics  of  the 
motomeuron.  Fluctuations  in  the  common  drive  are  reflected  in 
concurrent  fluctuations  in  the  firing  rate  of  motor  units  of  the 
same  pool. 

Furthermore,  it  has  been  shown  that  various  properties  of  a 
motor  unit  are  interrelated  [4,  6,  8].  For  instance,  Hcnneman’s 
well-known  size-principle  states  that  smaller  motomeurons  and 
equivalently  smaller  motor  units  become  active  at  lower  levels  of 
excitation  [8],  Smaller  motor  units  are  also  known  to  have  lower 
twitch  forces  with  longer  relaxation  times  than  larger  motor  units. 
Furthermore,  motor  units  with  lower  recruitment  thresholds, 
which  can  be  assumed  to  be  the  smaller  ones,  maintain  higher 
firing  rates  than  those  with  higher  recruitment  thresholds.  Yet 
another  relationship  exists  between  the  rank  or  recruitment 
threshold  of  the  motor  unit  and  its  firing  rate  response  to  increases 
in  excitation  [6],  In  summary,  the  recruitment  rank  of  the  motor 
unit  appears  to  define  many  of  its  firing  (electrical)  as  well  as 
twitch  (mechanical)  characteristics.  These  relationships  have  been 
exploited  in  defining  and  establishing  the  relationship  between 
various  characteristics  of  motor  units  in  the  present  model. 

B.  Firing  Rate  Determination 

Algorithms  to  produce  the  firing  rates  from  the  net  drive  for 
any  particular  MU  were  developed  based  on  the  force  firing  rate 
curves  reported  in  [6].  In  so  doing,  we  assumed  that  the  targeted 
force  level  represented  the  drive  to  the  muscle.  The  equations 
reported  in  [6]  along  with  curve  fitting  were  used  to  generate 
closed-form  representations  of  the  drive-firing  rate  characteristics 
for  each  motor  unit  based  on  its  recruitment  rank.  Figure  2 
represents  the  distribution  of  the  drive-firing  rate  profiles  among 
motor  units  of  the  pool. 


Integral  Pulse  Frequency  Modulation  (IPFM)  [2]  was  used 
to  generate  a  firing  train  with  a  firing  rate  equal  to  the  value 
determined  by  the  drive-firing  rate  profile  at  the  previous  stage  of 
the  model.  IPFM  is  essentially  an  integrate-and-fire  system  which 
integrates  the  input  over  time  until  a  threshold  (one  for  simplicity 
here)  is  reached.  At  this  point,  the  integrator  resets  itself  to  zero 
and  outputs  a  pulse  of  magnitude  one.  Hence,  a  lower  input  will 
result  in  firings  to  be  spaced  further  apart  compared  to  a  larger 
input.  The  operation  of  the  IPFM  model  is  exemplified  in  Figure 
3  for  two  different  input  levels. 
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Figure  3.  Generation  of  firing  trains  with 
Integral  Pulse  Frequency  Model. 


D.  Generation  of  Motor  Unit  Twitches 


In  determining  the  shape  of  the  force  twitch  the  motor  unit 
generates  in  response  to  a  firing  of  the  motomeuron,  Fuglevand,  et 
al.  [7]  developed  the  following  equations  after  [9]: 
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Figure  2.  Force  vs  firing  rate  relationships 
assumed  for  motor  units. 


I  lere,  f(l)  is  the  twitch  waveform  of  motor  unit  i.  It  depends  on 
the  parameters  Pt  and  Th  which  represent  the  twitch  peak  force 
and  contraction  time,  respectively.  The  coefficient  b  was  set  as 
b=(ln  RP)/n,  where  RP  is  the  range  of  peak  twitch  forces  and  n  is 
the  number  of  MUs.  TL  represents  the  longest  duration 
contraction  time  and  c  was  set  to  c=logRciRP ■  RCT  represents  the 
range  of  contraction  times  for  the  pool  of  MUs.  Contraction  time 
is  defined  as  the  time  from  zero  force  to  peak  force  [7],  These 
definitions  and  equations  were  adopted  without  change  in  the 
present  work.  Figure  4  displays  the  force  twitches  for  various 
motor  units  in  a  pool  of  150  motor  units.  In  this  simulation,  a 
maximum  contraction  time  (TL)  of  90  milliseconds,  a  contraction 
time  range  (RCT)  of  3,  and  a  peak  twitch  force  range  of  35  were 
assumed. 

In  defining  the  distribution  of  motor  unit  recruitment 
thresholds  we  again  adopted  the  relationships  described  by  [7]: 
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Figure  4.  Twitch  forces  for  motor  units  of 
various  ranks. 


RT (l)  —  ea  ‘  ,  where  RT  is  the  recruitment  threshold  and 
i  is  the  MU  number.  The  coefficient  a  establishes  a  range  of 
recruitment  thresholds.  It  was  set  to  a  =  (In  RR)/n,  where  RR  is 
the  maximum  recruitment  threshold  for  the  pool  and  n  is  the  total 
number  of  MUs  being  modeled  [7].  A  plot  of  recruitment 
threshold  versus  MU  number  is  shown  in  Figure  5  for  a  pool  of 
150  motor  units  with  a  recruitment  range  up  to  70  %MVC.  These 
values  were  selected  to  emulate  the  Tibialis  Anterior  muscle  on 
which  the  targeted  force  versus  firing  rate  values  were  based. 

III.  Validation  of  the  Model 

The  model  was  validated  using  three  methods.  First,  it  was 
shown  that  the  mean  firing  rates  of  MUs  from  the  model  matched 
the  characteristics  from  those  produced  by  real  MUs.  Secondly, 
the  power  spectrum  densities  of  compound  muscle  force  outputs 
from  the  model  were  compared  to  those  of  force  outputs  from 
human  muscle.  Thirdly,  the  output  of  the  model  when  driven  by  a 
ramp  input  was  verified.  Although  various  parameters  could  be 
fine-tuned  to  improve  the  results,  the  outputs  generated  by  the 
model  were  in  general  realistic.  Figure  6  represents  the  force 
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Figure  6.  Validation  of  the  model  based  on  the 
force  output  produced. 

output  generated  by  the  model  (with  parameters  as  defined  for 
Figure  4)  when  driven  with  a  linearly  increasing  input.  The  result 
approximates  a  linear  increase  that  would  be  expected  albeit  with 
some  deviation  in  the  middle  portion  of  the  excitation  range. 


IV.  SAMPLE  APPLICATIONS 


E.  Spike-Triggered  Averaging  to  Estimate  Twitches 

The  model  was  first  applied  to  the  validation  of  the  spike- 
triggered  averaging  (STA)  [3,  10]  technique  to  estimate  the  twitch 
force  of  individual  motor  units  from  the  compound  muscle  force 
output.  STA  estimates  the  force-time  characteristics  of  an 
individual  twitch  by  summing  and  averaging  multiple  time  frames 
aligned  with  the  individual  firings  of  the  motor  unit  in  question 
from  the  compound  muscle  force  output.  The  underlying 
assumptions  behind  STA  are  that  a  motor  unit’s  individual 
twitches  do  not  fuse  together,  and  the  twitches  of  other  motor 
units  are  uncorrelated  with  those  of  the  motor  unit  in  question. 
Fusion  of  twitches  would  hinder  the  contributions  other  than  the 
single  twitch  to  be  averaged  out  since  the  next  firing  of  the  same 
motor  unit  would  also  fall  in  the  averaged  frames.  Likewise,  the 
correlation  among  motor  unit  firing  patterns  would  keep  the 
contributions  of  other  units  from  falling  at  random  instances 
within  the  frames  taken  with  respect  to  the  firings  of  the  motor 
unit  under  investigation,  hence  prevent  them  from  canceling  out 
on  the  average. 

We  used  the  model  discussed  above  to  investigate  how 
confining  a)  higher  firing  rates  and  b)  correlation  among  the 
firings  of  motor  units  were.  Figure  7  presents  the  effect  of  firing 
rate  on  parameters  commonly  used  to  characterize  twitches:  peak 
force,  contraction  time  and  half-relaxation  time.  Model 
parameters  were  the  same  as  those  for  Figure  4.  The  motor  units 
were  driven  independently  (i.e.,  common  drive  was  zero)  in  order 
to  remove  the  effects  of  correlation  on  the  estimates.  It  can  be 
noted  that  even  at  firing  rates  that  are  too  low  to  be 
physiologically  viable,  the  estimation  error  reaches  unacceptable 
levels. 


Figure  5.  Distribution  of  motor  units  in  a  pool  of  The  effects  of  correlation  among  motor  unit  firing  patterns 

150  motor  units  with  recruitment  up  to  70  %  MVC.  were  investigated  by  varying  the  SNR  (defined  as  the  ratio  of  the 
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Figure  7.  Effects  of  motor  unit  firing  rate  on  of 
twitch  force  estimated  by  STA. 

signal  or  common  drive  to  noise).  Figure  8  presents  the 
estimation  results  for  the  three  twitch  parameters  at  SNR  values  of 
1,  2,  and  infinite.  The  model  parameters  were  the  same  as  before. 
The  drive  was  set  at  5%  MVC  in  order  to  minimize  the  effects  of 
fusion.  It  can  be  noted  that  even  at  SNR  values  of  1,  the  twitch 
estimation  becomes  so  noisy  that  the  derivation  of  twitch 
parameters  becomes  a  futile  effort. 
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Figure  8.  Effects  of  correlation  among  firing 
behavior  on  of  twitch  force  estimated  by  STA. 

F.  Hand  Dominance 

The  model  was  used  to  investigate  the  effects  of  hand- 
dominance  on  the  behavior  of  motor  unit  firing  rates.  In  specific, 
the  known  differences  among  the  dominant  and  nondominant 
muscles  were  implemented  and  simulations  were  run  to  determine 
what  could  be  concluded  about  the  firing  rates  or  drive  to  the  pool. 
Figure  9  presents  the  results  of  simulations  run  to  represent  the 
dominant  and  nondominant  muscles.  More  explicitly,  although 
slower  and  smaller  motor  units  were  assumed  for  the  dominant 
side  to  represent  the  shift  toward  more  Type  I  fibers,  their 
recruitment  range  was  narrower  to  parallel  experimental  findings 

[1],  The  number  of  motor  units  were  the  same  in  both  cases.  It 
can  be  seen  from  Figure  9  that  given  these  parameters,  for  the 
same  drive,  the  dominant  muscle  would  overshoot  the  target  if  it 
did  not  generate  lower  firing  rates  than  the  nondominant  side  as 
seen  in  experimental  data  [1],  It  remains  to  be  determined 
whether  the  lower  firing  rates  are  effected  via  a  lower  excitation 
or  changes  in  the  drive-firing  rate  characteristics. 

V.  CONCLUSIONS 

The  implemented  model  can  be  useful  in  studying  various 
techniques  or  physiological  paradigms.  It  goes  without  saying 
that  the  model  is  based  on  a  multitude  of  assumptions  and  the 
results  can  be  as  valid  as  the  assumptions  made.  In  general,  the 
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Figure  8.  Effects  of  correlation  among  firing 
behavior  on  of  twitch  force  estimated  by  STA. 

strength  of  modeling  studies  lies  in  not  proving  hypotheses  but 
disproving  unlikely  cases,  and  investigating  the  overall 
consequence  of  paradigms  with  conflicting  effects,  such  as  fatigue 
which  causes  a  decrease  in  twitch  amplitude  while  increasing 
twitch  duration.  The  availability  of  more  detailed  physiological 
information  will  improve  the  model,  as  will  the  inclusion  of 
factors  such  as  the  nonlinear  summation  of  twitches  and  time- 
dependence  in  the  model. 
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